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Abstract 


We propose a two-phase algorithm for solving continuous rank-one 
quadratic knapsack problems (R1QKPs). In particular, we study the so- 
lution structure of the problem without the knapsack constraint. In fact 
an O(nlogn) algorithm is suggested in this case. We then use the solu- 
tion structure to propose an O(n? logn) algorithm that finds an interval 
containing the optimal value of the Lagrangian dual of RIQKP. In the 
second phase, we solve the Lagrangian dual problem using a traditional 
single-variable optimization method. We perform a computational test 
on random instances and compare our algorithm with the general solver 
CPLEX. 
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1 Introduction 


The quadratic knapsack problem (QKP) deals with minimizing a quadratic 
function over one allocation constraint together with simple bounds on deci- 
sion variables. Formally, this problem can be written as 


minimize 


subject to @ 
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0<2<u4, (1c) 


where @ is a symmetric n x n matrix, a,¢,u € R”, and b € R. The 
QKP as a quadratic optimization problem is polynomially solvable when Q 
is positive definite matrix [12]. 

When Q is diagonal with strictly positive diagonal entries, then QKP can 
be viewed as a strictly convex separable optimization problem that has many 
applications (e.g., resource allocation [1, 13, 14] and multicommodity network 
flows [9]). The solution methods for solving this type of QKPs usually rely 
on the fact that the optimal solution to the Lagrangian dual subproblems 
can be explicitly obtained in terms of the Lagrange multiplier of (1b). 
Therefore, the problem reduces to find a value for such that the solution 
to the corresponding Lagrangian subproblem is satisfied equality constraint 
(1b). 

Helgason, Kennington, and Lall [9] proposed an O(n log n) algorithm for 
solving the equation based on searching breakpoints of the Lagrangian dual 
problem. Brucker [2] found an O(n) bisection algorithm based on the prop- 
erties of the Lagrangian dual function. Dai and Fletcher [4] proposed a two- 
phase method. A bracketing phase determines an interval containing the 
solution followed by the secant phase that approximates the solution within 
the promising interval. This method is modified by Comminetti, Mascaren- 
has, and Silva [3] by ignoring the bracketing phase and using a semi-smooth 
Newton method instead of the secant method. Liu and Liu [11] considered a 
special case of the strictly convex form of the problem. They found the solu- 
tion structure of the subproblems and used it in a modified secant algorithm. 

Robinson, Jiang, and Lerme [15] used the geometric interpretation of the 
problem and proposed an algorithm that works in the primal space rather 
than the dual space. This algorithm iteratively fixes variables and terminates 
after at most n iterations. 

In a more general case, when Q is positive semidefinite in (1), Dussault, 
Ferland, and Lemaire [7] proposed an iterative algorithm in which a QKP 
with diagonal Q should be solved in each iteration. Paradalos, Ye, and Han 
[12] suggested a potential reduction algorithm to solve this class of QKP. di 
Serafino et al. [6] proposed a two-phase gradient projection with acceptable 
numerical performance compared to similar gradient-based methods. 

QKPs with positive definite @ are also solved by a gradient projection 
method [4] and an augmented-Lagrangian approach [8}. 

In this paper, we suppose that Q is a rank-one symmetric matrix, that is, 
Q=qq' for some q € R”. Moreover, we assume that 0 < @. Without loss of 
generality, we assume that q; 4 0 for each 7. By the changing variables 

Li Ei, ci < am a; ¢ Lad uz <— max{0, gti; }, 
Gi Gi 


problem (1) is reduced to 
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1 
minimize 5 (Li a)” —e!'z, (2a) 
subject to ala =b, (2b) 
O<a<u (2c) 


Sharkey and Romeijn [16] studied a class of nonseparable nonlinear knap- 
sack problems in which one has to 


minimize g(s'x)—c'a, 
subject to a'a=b, (3) 
l<a<u, 


where g : R — R is an arbitrary real-valued function, and s € R” is given. 
They introduced an algorithm for solving (3) that runs in O(n? max{log n, ¢}), 
where ¢ is the time required to solve a single-variable optimization problem 
min{g(S)-—aS : L < $< U} for given a,L,U € R. With g(t) =?¢? and s 
equal to the all-one vector, problem (2) is a special case of problem (3). That 
is, there exists an O(n? max{logn,1}) = O(n? logn) algorithm for solving 
problem (2). 

In this paper, we consider a two-phase algorithm for solving problem (2). 
In Section 2, we study the solution structure of the relaxed version of the 
problem in which the equality constraint (2b) is excluded. We show that the 
relaxed version could be solved in O(nlogn) time. In Section 3, in phase I, 
we use the solution structure of the relaxed version to find an interval that 
may contain the optimal value of the Lagrangian dual function. This is done 
in O(n? logn) time in the worst case. Then, we perform phase II, in which 
we explore the interval by a single-variable optimization method to find the 
optimal Lagrangian multiplier with the desired precision. In Section 4, we 
perform a computational test. In particular, we compare the algorithm with 
the general quadratic programming solver CPLEX. 


2 Solution structure of the relaxed version 


In this section, we consider the following relaxed version of the problem (2): 


1 
minimize f(x) = 5 (lia) —c'a, (4a) 
subject to O<a<u. (4b) 


We propose a characterization of the solution in the primal space. Note 
that most of algorithms for such problems use the so-called KKT conditions 
to study the solution structure. 

Without loss of generality, we assume that cy > cg >-:: > cn = 0, 
and that 1; = 0,7 = 1,...,n. Given two vectors a,b € R”, we denote the 
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set {© : a <a < db} by [a,b]. Finally, given a vector u € R”, we define 
Up = yey Ui for k =1,...,n, and Up := 0. 
Now consider the following preliminary lemmas: 


Lemma 1. For k = 1,...,n define 1) as 


(k) Ug bee RS 
0, t=k+1,...,n, 


and x) as the all-zero vector, and define G;, as 
1 1 
Gr= 5 (Uk + Ug—1) — cp = Up_i + gk — Ck. 
Then the following assertions hold: 


(i) If % is the smallest index in {1,...,n} such that Gp > 0, then 


(ii) If Gy, <0 for all K=1,...,n—1, then minjey,. nn f(a) = f(a@™). 


Proof. (i) For 1 <k <n-—1, we have 


L 1 
Gy — Gry = 5 (Ux + Up_1) — Ce 5 (Un41 + Uk) + Choi 


1 
= 5 (up + Upsi) + (Ch+1 — CK) 


<0. 


Thus 
Gy < Go <-++ < Gp_1 <0 < Ga < Gag < +++ < Gp. 


On the other hand, for 1 < k <n—1, we have 


1 k+1 1 k 
(k+1)) _ (kK)) — =yy2 0 flee ef F2 ns 
f(« ) f@ ) — 5 Ui+1 2, UC 5 Uk + 2, Uj, 
1 1 
= sit — Uk41Ck41 — 5U 
1 
= 5 Vis = Uz) — Uk+1Ck+1 (5) 


I 


1 
9 (e+ — Ux) (Uni + Ue) — Un4ice41 


1 
= Uk41 (5 (Wis +Uz) — cxs1) 
= UnqiGe+1- 


Now let m >n—1. Then 
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PO ie Fa) Ee) ep ee) eee Tae) pe) 
= 0. 


Similarly, ifm </#—1, then f(a) — f(2@-Y) > 0. 
(ii) The second part can be easily proved by considering (5). 


We need the following result for two-dimensional version of problem (4). 


Lemma 2. Consider the following optimization problem: 


1 
minimize f(x1,x%2) = 5 (1 + £2)? — e121 — C22, 
subject to O0< a2, < wy, (6) 
0 < LQ < U2; 


where cy > cg > O and wu, and wz are real positive constants. Define set 
I := 1, Uh, where 1) = {(u1,42) : O< vq < ug}, and Ig = {(41,0) : O< 
x1 < wuz}. Then, problem (6) has no optimal solution outside of J. 


2 


2 + %2) = 32" — 
cz = g(z), where z = 41 + %. It is easy to see that «* = (xj, 23) with 
x} = min{c1,ui} and x} = min{c, — x}, uz}, is the optimal solution to the 
problem, and we have «* € I. Assume that c; 4 cg. The feasible region of (6) 
is equal to 1) UIn UI3 Uy, where [3 = {(21,22) : 0< a1 < u1,0 < rq < ug} 
and I, = {(0,42) : 0 < a < ug} U{(a1,u2) : O< a <u}. We show 
that there is no optimal solution in Jz; and I4. Indeed, we write the KKT 
optimality conditions as follows: 


Proof. If cy = cz, then f(a, 22) = 4(a#1 + 22)? — er (a1 


%1+%2-C1 +a, — a2 =), (7) 
Ly +%q— 2+ B — Bo =0, (8) 
ai(a1 — U1) =0, aga, = 0, (9) 
Pi(ag—U2) =0, Boxe =0, (10) 
O<a, <u, (11) 
0 < 22 < ua, (12) 
a, a2, G1, B2 = 0, (13) 


where q; and (6;, 1 = 1,2 are KKT multipliers corresponding to the bound 
constraints. If (#1, 22) € Is, then from (9) and (10), we have ay = a2 = $y = 
(G2 = 0. Substituting these values in (7) and (8) implies that c, = c2, which 
contradicts our assumption. On the other hand, if (#1,2%2) € I4 and x; = 0, 
then a, = 0. Now, (7) implies that 2 = cy + a2 > 0. Thus 62 = 0. From 
(8), we have wg = co— (31. Therefore, cp = cy +a2+1 > cy. This contradicts 
our assumption on ¢;’s. That is, problem (6) has no optimal solution with 
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x; = 0. Now, if v2 = ue, then ay = 0 and £2. It implies from (7) and (8) 
that 0 < ag + 6; =cg—c, <0. That is, co = c,, a contradiction. 


Theorem 1. Suppose that 2) and G;, k =1,...,n, and fi are defined as 
in Lemma 1. Then the following assertions hold: 


(i) For 7 > 1, define 6; and 62 as 6, := min {cp_1 — Un—2, un_i} and 
d2 := max {cq — Un_1,0}. Also, define % and & as 


where e; is the ith column of the identity matrix of dimension n. Then 
min{ f(Z), f(Z)} is the optimal value of the following optimization prob- 


lem: ee 
minimize x), 
fle) . - 
subject to 2) <a <a), 


(ii) For 7 = 1, define 6 := min{c),ui} and & := de;. Then f(%Z) is the 
optimal value of the following optimization problem: 


minimize f(x), 


subject to 2 <a <a. 


(iii) For G, < 0 for all k = 1,...,n, define 6 := min{c, — Un—1, un} and 
& := a") + be,. Then f(%) is the optimal value of the following 
optimization problem: 


minimize f(a), 


subject to c"Vcg<g™, 


Proof. (i) By Lemma 2, we can partition the optimal solution set as I; U Jo, 
where I, = [x~?), 2-0] and Ip = [2-2]. We show that f(z) = 
min{ f(x) : « € Ly} and f(z) = min{f(ax) : « © Ig}. Indeed, we use a 
simple technique of single-variable calculus. Let « € I). Then x = 2(6), 
for some 6 € [0,un—1], where 2(5) = x'"-?) + den_1. Thus the problem 
min{ f(z) : x € I;} reduces to min{ f(x(d)) : 0 < 6 < upj_i}. On the other 
hand, one can write 


n-2 
f(x(6)) = ; (Un—2 + 6) _ > Cj; — Ch—10. 
i=l 


We have df(x(0))/dd = Un—2 + 6 — cn_1. Thus df(x(6))/dd = 0 only if 
6 = 6' =ca_1 — Up—z. Since d? f(x(5))/dd? > 0 and 6’ > Suq_1 the optimal 
value is achieved at 61. 

To prove f(%) = min{f(x) : x € Iz}, by the same argument as 
the previous paragraph, it suffices to solve single optimization problem 
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min{ f(a(5)) : 0 < 6 < un}, where 2(5) = a-) + den. It is easy to see 

that if 6 = 6 = cg, — Up_1, then df(a(5))/dd = 0. Since 6’ < Sua, by the 

definition of 7, then f(%) is the optimal value of min{f(«) : « € Ip}. 
The proof of parts (ii) and (iii) is similar. 


The following Corollary 1 presents simple conditions under which the 
optimal solution to the problem in Theorem 1(i) is or @. 


Corollary 1. In Theorem 1(i), 
(i) if 6) = unj_i, then min{ f(Z), f(%)} = f(@), and 
(ii) if dg = 0, then min{ f(z), f(z)} = f(z). 


Proof. For brevity, we just prove part (i). The proof of the second part is 
similar. Under the assumption of part (i), we have 


n—-1 


f(z)- f(z) = UF 1 > CyUi— 5c at), Cus ten (c Cah —Un-1) = ; (Un-1 = cn) > 0. 


w=1 


Theorem 1 solves a relaxed version of problem (2). In Theorem 2, we 
show that the solution to the relaxed version is the solution to the original 
problem. 


Theorem 2. Define G;’s, x“)’s, 7, %, and % as in Theorem 1. Then, the 
following assertions hold: 


(i) If1<nv<n, then min{f(Z), f(Z)} is the optimal value of (4), where @ 
and & are defined as in Theorem 1(i). 


(ii) If mi =1, then f(Z) is the optimal value of (4), where & is defined as in 
Theorem 1 (ii). 


(iii) If Gy < 0 for all k =1,...,n, then f(Z) is the optimal solution to (4), 
where ¢ = «(—)) 4+ 6’e, and 6’ = minfe, — Up_1, Un}. 


Proof. For two vectors x, z € R”, we have 


fle) — fle) = 5172 +17 aa 1'x)—c'(z-2). (15) 


Let x be a feasible solution to (4). If 2 = u =a, then from the defini- 
tion of n, we have f(x"~!) < f(x), and the result follows from Theorem 1. 
Suppose « 4 u. We show there exists a specially structured feasible solution 
x’ that is better than x. Indeed, let k be such that 


U;, < ile < Uxg41- 
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Define vector x’ by 


Ui; = ey oe 
ge=<1'x-U,, i=k+1, 
0, t=k+2,...,n. 


Then, clearly wx’ is feasible for (4) and 


k n k n k 
Wek ! / Ss oh * o > T 
a=) UTC «;=> wi + > «i — uj=1 cz. 
i=l i=1 i=l i=1 


t=k+2 


Moreover, we obtain 


k k n k 
cla’ = S- Use, + Ch41 0p 41 = S- UWiCis + Ck41 S- Xi — Ck4+1 S- Uj 
i= i=1 i=1 i=1 
k k n k 
= S- UiCi + Ck+1 x Ui + Ce+1 S- Li — Ck+1 pa U4 
i= i=l i=k+1 i=1 
k k n 
> S- WiC; + So (ai —u)o t+ S- x;c; (by the monotonicity of c;’s) 
i= i=l i=k+1 
=c'z 


Therefore, (15) implies that f(x’) — f(x) = —c' (x' — 2x) < 0, that is, f(z’) < 
f (2). 

(i) Now, we consider three cases for the index & introduced in the defini- 
tion of a’: k>n,k<n—-2,andk=n—1,n—2. In the latter case, we have 
a2) < a! < «™), so the assertion is true by Theorem 1, since 


min{ f(@), f(@} = minf f(x) : «€ fe), 2™)} < f(a’) < f(a). 


better than 2’, that is, f(2) < f(x’) for some i = 1,...,n. By Lemma 1, 
f(@@-Y) < f(a) and the result follows by Theorem 1. 


First, let k > n. Then 


eR 


f(a) _ f(z’) = 7A To (k) _ 1 e'\(1' 2 a 1'2') = cl (a) = a’) 


1 
= ~ 5 te+1 2Ue + Ue41) a5 Cho 12h 41 


1 
Tet (S0n + Chan) cna): 


I 


On the other hand, we have 
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n-1 k 
QU + Lies =230 ut > ut ee > e+ Ue. 
i=1 i=n 
Therefore, 
1 F 1 
5 (2Ui t p41) Ck+1 oe 3 (Un + Un-1) — cn = Ga 2 0. 


Thus f(x) < f(z’). 
Now, let k <n — 2. Then 


f(a) — f(a’) = 


(aT aD) Te (aT et) 417 2!) — eh (et) — 2’) 


1 
2 
1 
5 (Uet1 — Tegr)(Uner + Ur + tir) — Ch+i (tert — Tey) 


1 
= (Uk+1 — Th41) (S0u + B41 + Uti) — cxs1) 


On the other hand, we have 


n—2 
2Up 4+ 241 tnt < 2U_ 4+ ued. < 2U, +2 S- Uy +Un—1 = Un—2 + Up-1. 
i=k+1 


Hence, 


1 
5 (Un—-2 + Un ) Ca 1=Gr 1 <0. 


L / 
5 (2Uk t Terr t Uk+1) Ck+1 < 


That is, f(x*t)) < f(x’). Thus in both cases, there exist a point, say x, 
such that f(x) < f(a’) < f(a). Now, by Lemma 1, f(a#(®-)) < f(@) < 
f(x) and the result follows by Theorem 1. 

Proof of (ii). Consider the possible values of k at the beginning of the 
proof of part (i). Here, we just have k > m = 1. Now, similar argument for 
this case proves (ii). 

Proof of (iii). Again consider the possible values of & at the beginning of 
the proof of part (i). Similar argument with case k < i —2 form =n+1 
proves part (iii). 


We conclude the following result on the time needed to solve problem (4). 
Theorem 3. There exists an O(n logn) time algorithm for problem (4). 


Proof. When the index 7 is determined, the solution can be determined in 
O(n) time. We need O(nlogn) to sort the vector c, O(n) to compute vector 
G, and O(log n) to find the index i. That is, problem (4) can be solved in 
O(nlogn). 
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3 The algorithm 


In this section, we propose our algorithm for solving problem (2). The al- 
gorithm consists of two phases: bounding the optimal Lagrangian multiplier 
and computing the optimal solution to the desired precision. The bounding 
phase is based on the Lagrangian dual of (2) and the solution structure of 
the relaxed version has been described in section 2. 


3.1 Lagrangian dual 


Let A be the Lagrange multiplier of equality constraint in (2). Then, the 
Lagrangian function is given by 


P(A) = min 5(7 ey? —c'a+X(b—-a'z) : 0<a< ul 
(16) 
= do -nin 517) —(c+)a)"s > 0<a2< ul : 


We have the following statement about the structure of the Lagrangian 
function ¢. 


Theorem 4. For a given Lagrange multiplier A, define 7 as in Theorem 2. 
If ni > 1, then 


G(A) = ABt fa(w@O), if en(A) S Un-1 S en-1A), (Type A) 
(A) = Ab + pa(A), if Un—1 < en(A), (Type B) 
P(A) = AB + pa—1(A), if Un—1 > cn—1 A), (Type C) 


where fy is the objective function of the optimization part of (16), and 


1 1 
pR(A) = = 50” ~@) aed 5 ck ees 
dy = glk) + (cK _ Up_-1)€k- 


Proof. The proof is based on the four possible cases for 6; and 62 in The- 
orem 2. We just prove (Type A) and, for the sake of brevity, we omit the 
remaining parts. 

Suppose that ca(A) < Ua—1 < ca_i(A). Then we have cq_1(A) — Un—2 > 
Un—1 and ca(A) — U;_1 < 0. Therefore, the values of 6; and 62 in Theorem 2 
can be determined as 


61 =min{ca—1(A) — Un—2, Ua—-1} = Ua-1, 
d2 = max{ca (A) = ORa33 0} = ca(A) a Up_-1.- 
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Thus we have @ = 2("—-2) + dye,_1 = cl), By some simplifications, we 
have fy(#) — f,(Z) = $(cn(A) — Un_i)? > 0. Now, Theorem 2 implies that 


min{ f(z) : 0< a <u} = min{ f(z), f(Z)} = fa(®) = fala). 


Now, one may conclude that if 7% is fixed on an interval [Ag, A»], then 
@(A) is a piecewise function that contains exactly three pieces. However, the 
following simple example shows that this is not true. 


Example 1. Consider problem (2) with the following parameters: 


a’ =[-7-57-57], c' = [54 4415-8 —70], 
ul = [62 48 36 84 59]. 


In Figure 1, we plot d(A) for A € [—8.36, 7.00]. We distinct three cases in 
(Type A), (Type B), and (Type C) in blue, red, and green, respectively. As 
it can be seen in Figure 1, @(A) consists of four pieces. 


-1000 


-1500 - 


-2000 


S} -2500 


-3000 > 


-3500 - 


@ values of type | 
@ values of type II 
¢ values of type III 


-10 -8 -6 -4 -2 0 2 4 6 8 
—8.36 < A < 7.00 


Figure 1: Plot of the Lagrangian function ¢(A) for Example 1. 


The inner optimization problem in (16) is a special case of problem (4) 
that can be solved by Theorem 2. In Theorem 2, it is assumed that coefficients 
of the linear term in the objective function are sorted in decreasing order. In 
problem (16), the order of coefficients of the linear term depends on ». From 
now on, we denote by c;(A) the coefficient of x;, that is, c;(A) = cj + Aaj. 
Moreover, we denote the line {c;(A) : A € R} by 2;. It is easy to see that 
when 2 becomes greater than the intersection of ¢; and £;, c;(A) and c¢;(A) 
change position in the ordered list of coefficients. 

We use a modification of the well-known plane sweep algorithm to find the 
ordered intersection points of lines {€; : i=1,...,n}. Now, let Aq, ’, and 
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A» be three consecutive intersection points. Then, because the Lagrangian 
function is unimodal, the optimal Lagrange multiplier \* lies in the interval 
[Aa As] if O(A’) > O(Aa) and 6(0’) > oro). 

We modify the implementation of the line-sweep algorithm proposed in 
[5]. In this algorithm, a vertical line @ sweeps the plane from left to right. 
The status of the sweep line is the ordered sequence of lines that intersect it. 
The status initially contains all lines in the order of decreasing slope, that is, 
the order of lines when they intersect with the sweep line at 1 = —oo. The 
status is updated when ¢ reaches an intersecting point. For example, suppose 
that the sequence of four lines &, ¢;, €;, and €, appears in the status when ¢ 
reaches the intersection point of ¢; and @;. Then, £; and €; switch the position 
and intersection of lines £; and @,,, and the intersection of 0; and & are to 
be checked. The new detected intersection points are stored to proceed. The 
order of cost coefficient of the linear term in ¢(A) is unchanged between two 
consecutive intersection points. 

If c;(A) < 0 for some 7, then x; = 0 in the optimal solution to the ¢(A) 
subproblem. We introduce a set Z to store the non-vanished variables. To do 
so, we add a dummy line 0g : cg(A) = 0. In each intersection of the dummy 
line and the other lines, the set Z should be updated. In fact, if 2; intersect £9 
andi ¢ Z, then we add i to Z; otherwise, if 7 € Z, then it should be removed 
from Z. In other words, since we sweep the plane from left to right, if @; 
intersect &p and a; < 0, then we addi to Z. If @; intersect @p and a; > 0, then 
it means that 7 should be removed from Z. Moreover, Z initially contains the 
set of all lines with a positive slope. With this modification, we guarantee that 
between two consecutive intersection points, the set of zero-valued variables 
is unchanged. It should be noted here that lines with equal slopes are sorted 
based on increasing order of c;’s. We summarize the approach in Algorithm 
1. This algorithm is used as the first phase in the main algorithm. 


Theorem 5. Algorithm 1 runs in O(n? log n) time. 


Proof. Initializing state array , line indices array p and the queue Q in 
steps 3-7 needs O(nlogn) time. In each iteration, we perform two main 
operations: computing the value of ¢ for a new intersect point A"°Y and 
updating Q. The order of c;(A"®”) and the vector G can be updated from 
the previous intersection point in O(1) time. Finding 7 needs O(log n), using 
binary search. On the other hand, insertion and deletion on the priority queue 
Q takes O(log) since one can implement the priority queue by a heap to 
store the intersection points. Therefore, each iteration of the main loop needs 
O(log n) time. Since there are O(n?) intersection points, the algorithm runs 
in O(n? log n). 


Let Arg and Avg be the smallest and greatest intersection points of lines 
{é@; : «= 1,...,n}, respectively. The optimal solution to the Lagrangian 
problem may lie out of the interval [Arp,Auvp]. In this case, Algorithm 1 
fails to find the optimal interval. So, we explore the outside of [Arg, Aug] in 
a separate phase. 
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Algorithm 1 A plane sweep algorithm for finding an interval containing the 
optimal solution to the Lagrangian dual problem. 


a 
: Output: interval [\,, A] that contains the optimal solution to the prob- 


Input: vectors c, a, and u and scaler b. 


lem max,er (A) or the smallest and largest intersection points. 
Initialize a state array, @ = [1,...,n] , with lines ¢[1],..., [n] sorted in 
decreasing order of their slope. 
Initialize queue Q = @. 
Initialize line indices array p = [1,...,n]. 
FAIL ¢+ true 
Insert intersection points of all adjacent lines into Q. 
Set \prev — —0O, Prev prev = —oo 
while Q is not empty do 
Pop from Q the current intersection point 
two adjacent lines ¢[7] and ¢[,]. 
Update state array: £[p[t]] © €[p[y]]. 
Update the line indices array: pli] © ply]. 
Insert the intersection point of ¢[p[i]] and ¢[p[¢]+ 1] and the intersection 
point of é[p|j]] and £[p[j] — 1] into Q, if there exists any. 
if (APrev) > p(APrev Preys) and (APrev) > o(Are”) then 
Set FAIL «+ false 
return [Prev Prev \rew) as the promising interval. 
end if 
Set \prev prev oa Prev | 
Set APrev < Anew, 
end while 


Av©~ and the corresponding 


: if Fart then 


return the smallest A; and the largest Ayg intersection points. 
end if 
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First, consider the exploration of (—oo, Az). Since the components of 
vector G in Theorem 2 are linear functions in term of A, then there exists 
Nig < AxB such that the order of G;,’s does not change for \ < \/, 7. Indeed, 
a similar procedure for finding the smallest intersection of lines @;’s can be 
used here to compute A), 2. Now, since ¢(A) is unimodal, one can conclude 
that 

max (A) 


(—00,ALB] IAL p ALB] 


(A). (17) 


Similarly, for the values of A > Avg, one can find a threshold, say Aj, p, 
such that 
max @(A)= max (A). (18) 


[AuB,co) us Aye) 


We summarize the main algorithm in Algorithm 2. 


Algorithm 2 A two-phase algorithm for solving rank-one QKP (2). 


1: Run Algorithm 1 to find a promising interval that contains the optimal 
Lagrange multiplier. 

2: if Algorithm 1 returns an interval [\,,»] then 

Solve the optimization problem max ),,) (A) and return the optimal 

solution. 

: else 

Solve optimization problems (17) and (18) and store the optimal values. 

: end if 

: return the best \ found as an optimal Lagrange multiplier. 


It is clear that Algorithm 2 converges to the optimal solution since the 
output interval of Algorithm 1 contains the optimal solution and ¢()) is 
unimodal. In fact, the single variable optimization problem in step 3 can be 
solved efficiently by a classical root-finding algorithm. 


4 Computational experiments 


In this section, we compare the running time of Algorithm 2 with the general 
convex quadratic programming solver, CPLEX. We implement Algorithm 2 
with MATLAB R2019b. All runs are performed on a system with a Core 
i7 2.00 GHz CPU and 8.00 GB of RAM equipped with a 64bit Windows 
operating system. We solve single variables optimization problems (17), (18), 
and step 3 in Algorithm 2, using MATLAB built-in function fminbnd, which 
is based on the golden section search and parabolic interpolation. 

Our testbed contains two types of randomly generated rank-one knapsack 
problems up to n = 100,000 variables. In the first type, the vectors a and c 
are integral and generated uniformly from the same interval. We denote this 
type by Typel. In the second type (Typell), the vectors a and c are positive 


IJNAO, Vol. 12, No. 3 (Special Issue), 2022, pp 567-584 


A two-phase method for solving continuous rank-one quadratic ... 581 
Table 1: Parameters for two types of problem instances. 


Type a Cc l u—_l 


Typel U(—50,50)  U(—50,50) U(0,20) U(1,100) 
Typell U(—100,10) U(10,100) U(0,20) U(1,100) 


and negative randomly generated integral vectors, respectively. In Table 1, 
we summarize the parameter values for problem instances. 

As a well-known general convex quadratic programming solver, we chose 
CPLEX (ver. 12.9) to compare our results. 

Based on our numerical results, we set the quadratic programming solver 
of CPLEX (qpmethod option) to the barrier. The barrier convergence toler- 
ance, convergetol, is set to le— 12 (The default value is convergetol = le —6). 
It should be noted here that this setting is applied after we found that the 
default value leads to the optimal solutions that have components with a 
“meaningful” distance to their correct values. Another point is that for other 
optimizers such as primal and dual, CPLEX found the optimal solution in 
correct precision, but the running time is too long for large instances. For 
brevity, we do not report details related to these experiments. 

After completing our experimental tests, we found in [10] that the sparsity 
of the Hessian matrix influences the performance of CPLEX. To increase the 
performance, we reformulate our problem as 


1 
min { 5a? — 0% : 1's-—y=0, a's=b, o<ncul. 


We denote the results corresponding to running CPLEX on the original prob- 
lem and the aforementioned modification, respectively, by CPLEX,,, and 
CPLEX ef. 

Table 2 shows the average running time for five runs of each algorithm/- 
solver for each problem size. Dash sign, “-”, denoted the algorithm/solver 
encounters out-of-memory status. 

In all cases, Algorithm 2 outperforms CPLEX>,;g. For instances up to 
n = 5000, our algorithm and CPLEX, have the same running time, whereas 
for larger instances, CPLEX,¢f has smaller running time. 


5 Conclusions 


In this paper, we proposed a two-phase algorithm for rank-one QKPs. To this 
end, we studied the solution structure of the problem when it has no resource 
constraint. Indeed, we proposed an O(nlogn) algorithm to solve this prob- 
lem. We then used the solution structure to propose an O(n? log n) line-sweep 
algorithm that finds an interval that contains the optimal Lagrange multi- 
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Table 2: A comparison of running times (in seconds) between our algorithm 
and CPLEX yp, and CPLEX ye. 


n Our algorithm CPLEX,,g CPLEX;ef 
1000 Typel 0.06 0.09 0.01 
Typell 0.01 0.06 0.02 
1500 Typel 0.04 0.15 0.02 
Typell 0.02 0.13 0.02 
2000 Typel 0.04 0.27 0.02 
Typell 0.02 0.27 0.02 
5000 Typel 0.09 2.21 0.02 
Typell 0.06 2.12 0.03 
10000 Typel 0.26 16.26 0.04 
Typell 0.23 16.95 0.05 
15000 Typel 0.62 61.20 0.10 
Typell 0.63 65.88 0.10 
20000 Typel 1.16 - 1.20 
Typell 0.88 - 1.02 
50000 Typel 3.22 - 0.11 
Typell 3.19 - 0.11 
100000 ~—Typel 12.19 - 0.14 
Typell 11.31 - 0.17 
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plier. Then, the estimated optimal interval was explored for computing the 
optimal solution with the desired accuracy. Our computational tests showed 
that our algorithm has better running time than CPLEX when CPLEX is 
used to solve the original problem. After a reformulation of the problem, 
CPLEX outperforms our algorithm for instances with n > 5000. 
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